Companion_H24_5aDay
Clémentine Cottineau-Mugadza, Julien Perret, Romain Reuillon, Sébastien Rey-Coyrehourcq, Julie Vallée
2024-09-12
Libraries
library(tidyverse)
library(reshape2)
library(sf)
library(spdep)
library(tmap)
#library(OasisR) ## old package (OK for version R.2.1)Synthetic population - H24 library (section 2.1)
Import data for synthetic population located in ‘night’ cell
pop0, in ‘day’ cell pop1 and in ‘evening’ cell
pop2
pop_0 <- st_read( "data/population.gpkg",
layer="Population_0", quiet = TRUE)
pop_1 <- st_read("data/population.gpkg",
layer="Population_1", quiet = TRUE)
pop_2 <- st_read("data/population.gpkg",
layer="Population_2", quiet = TRUE)Create new colums : gender (2 groups) ; Age (3 groups); Educ (3 Groups)
pop_0 <- pop_0 %>%
rowwise() %>%
mutate(
men = sum(across(starts_with('pop_1'))),
women = sum(across(starts_with('pop_2'))),
pop = men + women,
age1 = sum(across(starts_with(c('pop_1_1', 'pop_2_1')))),
age2 = sum(across(starts_with(c('pop_1_2', 'pop_2_2')))),
age3 = sum(across(starts_with(c('pop_1_3', 'pop_2_3')))),
educ1 = sum(across(ends_with('_1'))),
educ2 = sum(across(ends_with('_2'))),
educ3 = sum(across(ends_with('_3'))),
pct_educ3 = (educ3 / pop) * 100,
pct_age3 = (age3 / pop) * 100,
pct_women = (women / pop) * 100
)
pop_1 <- pop_1 %>%
rowwise() %>%
mutate(
men = sum(across(starts_with('pop_1'))),
women = sum(across(starts_with('pop_2'))),
age1 = sum(across(starts_with(c('pop_1_1', 'pop_2_1')))),
age2 = sum(across(starts_with(c('pop_1_2', 'pop_2_2')))),
age3 = sum(across(starts_with(c('pop_1_3', 'pop_2_3')))),
educ1 = sum(across(ends_with('_1'))),
educ2 = sum(across(ends_with('_2'))),
educ3 = sum(across(ends_with('_3')))
)
pop_2 <- pop_2 %>%
rowwise() %>%
mutate(
men = sum(across(starts_with('pop_1'))),
women = sum(across(starts_with('pop_2'))),
age1 = sum(across(starts_with(c('pop_1_1', 'pop_2_1')))),
age2 = sum(across(starts_with(c('pop_1_2', 'pop_2_2')))),
age3 = sum(across(starts_with(c('pop_1_3', 'pop_2_3')))),
educ1 = sum(across(ends_with('_1'))),
educ2 = sum(across(ends_with('_2'))),
educ3 = sum(across(ends_with('_3')))
)Distribution of the synthetic population by demographic category in the Paris Region (night cells) (Figure 1)
Total Population density
bks_pop <- c(1, 10, 100, 1000, 10000, Inf)
tm_shape(pop_0) + tm_fill(col="pop", style="fixed",
breaks = bks_pop, palette="Greys") +
tm_legend(outside=TRUE)
Proportion of people with 2+ years of higher education
bks_edu <- c(0, 3.5, 7.7, 11.4, 17.2, 100)
tm_shape(pop_0) + tm_fill(col="pct_educ3", style="fixed",
breaks = bks_edu, palette="Reds") +
tm_legend(outside=TRUE)
Proportion of people aged 60+
bks_age <- c(0, 12.5, 17.1, 20.6, 25.1, 100)
tm_shape(pop_0) + tm_fill(col="pct_age3", style="fixed",
breaks = bks_age, palette="Blues") +
tm_legend(outside=TRUE)
Proportion of women
bks_sex <- c(0, 45.8, 50, 51.6, 54.5, 100)
tm_shape(pop_0) + tm_fill(col="pct_women", style="fixed",
breaks = bks_sex, palette="Greens") +
tm_legend(outside=TRUE)
Duncan Dissimilarity Index (Figure 1 and Appendix 3)
Duncan’s segregation index is one-group form of dissimilarity index DI and measures the unevenness of a group distribution compared to the rest of the population. It can be interpreted as the share of the group that would have to move to achieve an even distribution compared to the rest of the population. You should not include a column with total population, because this will be interpreted as a group.
In night cells
dfpop_0 <- pop_0 %>% st_drop_geometry()
pop_0_Duncan_18categ<- as.data.frame(ISDuncan(dfpop_0[,2:19]))
colnames(pop_0_Duncan_18categ)[1] <- "Duncan_night"
rownames(pop_0_Duncan_18categ) <- paste (colnames(dfpop_0[ , 2:19]))
pop_0_Duncan_sex<- as.data.frame(ISDuncan(dfpop_0[,20:21]))
colnames(pop_0_Duncan_sex)[1] <- "Duncan_night"
rownames(pop_0_Duncan_sex) <- paste (colnames(dfpop_0[ , 20:21]))
pop_0_Duncan_age<- as.data.frame(ISDuncan(dfpop_0[,22:24]))
colnames(pop_0_Duncan_age)[1] <- "Duncan_night"
rownames(pop_0_Duncan_age) <- paste (colnames(dfpop_0[ , 22:24]))
pop_0_Duncan_educ<- as.data.frame(ISDuncan(dfpop_0[,25:27]))
colnames(pop_0_Duncan_educ)[1] <- "Duncan_night"
rownames(pop_0_Duncan_educ) <- paste (colnames(dfpop_0[ , 25:27]))
pop_0_all_Duncan <- bind_rows(pop_0_Duncan_18categ, pop_0_Duncan_sex, pop_0_Duncan_age, pop_0_Duncan_educ)
pop_0_all_Duncan_round <- pop_0_all_Duncan %>% mutate(across(starts_with("Duncan"), round, 3))In day cells
dfpop_1 <- pop_1 %>% st_drop_geometry()
pop_1_Duncan_18categ<- as.data.frame(ISDuncan(dfpop_1[,2:19]))
colnames(pop_1_Duncan_18categ)[1] <- "Duncan_day"
rownames(pop_1_Duncan_18categ) <- paste (colnames(dfpop_1[ , 2:19]))
pop_1_Duncan_sex<- as.data.frame(ISDuncan(dfpop_1[,20:21]))
colnames(pop_1_Duncan_sex)[1] <- "Duncan_day"
rownames(pop_1_Duncan_sex) <- paste (colnames(dfpop_1[ , 20:21]))
pop_1_Duncan_age<- as.data.frame(ISDuncan(dfpop_1[,22:24]))
colnames(pop_1_Duncan_age)[1] <- "Duncan_day"
rownames(pop_1_Duncan_age) <- paste (colnames(dfpop_1[ , 22:24]))
pop_1_Duncan_educ<- as.data.frame(ISDuncan(dfpop_1[,25:27]))
colnames(pop_1_Duncan_educ)[1] <- "Duncan_day"
rownames(pop_1_Duncan_educ) <- paste (colnames(dfpop_1[ , 25:27]))
pop_1_all_Duncan <- bind_rows(pop_1_Duncan_18categ, pop_1_Duncan_sex, pop_1_Duncan_age, pop_1_Duncan_educ)
pop_1_all_Duncan_round <- pop_1_all_Duncan %>% mutate(across(starts_with("Duncan"), round, 3))In evening cells
dfpop_2 <- pop_2 %>% st_drop_geometry()
pop_2_Duncan_18categ<- as.data.frame(ISDuncan(dfpop_2[,2:19]))
colnames(pop_2_Duncan_18categ)[1] <- "Duncan_evening"
rownames(pop_2_Duncan_18categ) <- paste (colnames(dfpop_2[ , 2:19]))
pop_2_Duncan_sex<- as.data.frame(ISDuncan(dfpop_2[,20:21]))
colnames(pop_2_Duncan_sex)[1] <- "Duncan_evening"
rownames(pop_2_Duncan_sex) <- paste (colnames(dfpop_2[ , 20:21]))
pop_2_Duncan_age<- as.data.frame(ISDuncan(dfpop_2[,22:24]))
colnames(pop_2_Duncan_age)[1] <- "Duncan_evening"
rownames(pop_2_Duncan_age) <- paste (colnames(dfpop_2[ , 22:24]))
pop_2_Duncan_educ<- as.data.frame(ISDuncan(dfpop_2[,25:27]))
colnames(pop_2_Duncan_educ)[1] <- "Duncan_evening"
rownames(pop_2_Duncan_educ) <- paste (colnames(dfpop_2[ , 25:27]))
pop_2_all_Duncan <- bind_rows(pop_2_Duncan_18categ, pop_2_Duncan_sex, pop_2_Duncan_age, pop_2_Duncan_educ)
pop_2_all_Duncan_round <- pop_2_all_Duncan %>% mutate(across(starts_with("Duncan"), round, 3))Build table
pop_all_Duncan_round <- bind_cols(pop_0_all_Duncan_round, pop_1_all_Duncan_round, pop_2_all_Duncan_round)
pop_all_Duncan_round Moran’s index of spatial autocorrelation (Figure 1 and Appendix 3)
In night cells
dfpop_0_prop <- dfpop_0/dfpop_0$pop
pop_0_geom <- pop_0[, -1:-28]
pop_0_prop <- bind_cols(pop_0_geom, dfpop_0_prop)
nb0 <- poly2nb(pop_0_prop, queen=TRUE)
lw0 <- nb2listw(nb0, style="W", zero.policy=TRUE)
list <- colnames(dfpop_0_prop[,2:27])
pop_0_moran <- data.frame(matrix(ncol = 2))
colnames(pop_0_moran) <- c("pop_night_IMoran" , "pop_night_Moranpvalue")
for (i in list)
{
pop_0_i <- pop_0_prop[,i]
colnames(pop_0_i)[1] ="ok"
pop_0_moran[i, 1] <- moran.test(pop_0_i$ok, lw0, zero.policy=TRUE, alternative="greater") [3]
pop_0_moran[i, 2] <- moran.test(pop_0_i$ok, lw0, zero.policy=TRUE, alternative="greater") [2]
}
pop_0_moran <- pop_0_moran[-1,] # Delete empty first line
pop_0_moran_round <- pop_0_moran %>% mutate(across(starts_with("pop"), round, 3))
pop_0_moran_roundIn day cells
dfpop_1_prop <- dfpop_1/dfpop_1$pop
pop_1_geom <- pop_1[, -1:-28]
pop_1_prop <- bind_cols(pop_1_geom, dfpop_1_prop)
nb1 <- poly2nb(pop_1_prop, queen=TRUE)
lw1 <- nb2listw(nb1, style="W", zero.policy=TRUE)
list <- colnames(dfpop_1_prop[,2:27])
pop_1_moran <- data.frame(matrix(ncol = 2))
colnames(pop_1_moran) <- c( "pop_day_IMoran" , "pop_day_Moranpvalue")
for (i in list)
{
pop_1_i <- pop_1_prop[,i]
colnames(pop_1_i)[1] ="ok"
pop_1_moran[i, 1] <- moran.test(pop_1_i$ok, lw1, zero.policy=TRUE, alternative="greater") [3]
pop_1_moran[i, 2] <- moran.test(pop_1_i$ok, lw1, zero.policy=TRUE, alternative="greater") [2]
}
pop_1_moran <- pop_1_moran[-1,] # Delete empty first line
pop_1_moran_round <- pop_1_moran %>% mutate(across(starts_with("pop"), round, 3))
pop_1_moran_roundIn evening cells
dfpop_2_prop <- dfpop_2/dfpop_2$pop
pop_2_geom <- pop_2[, -1:-28]
pop_2_prop <- bind_cols(pop_2_geom, dfpop_2_prop)
nb2 <- poly2nb(pop_2_prop, queen=TRUE)
lw2 <- nb2listw(nb2, style="W", zero.policy=TRUE)
list <- colnames(dfpop_2_prop[,2:27])
pop_2_moran <- data.frame(matrix(ncol = 2))
colnames(pop_2_moran) <- c( "pop_evening_IMoran" , "pop_evening_Moranpvalue")
for (i in list)
{
pop_2_i <- pop_2_prop[,i]
colnames(pop_2_i)[1] ="ok"
pop_2_moran[i, 1] <- moran.test(pop_2_i$ok, lw2, zero.policy=TRUE, alternative="greater") [3]
pop_2_moran[i, 2] <- moran.test(pop_2_i$ok, lw2, zero.policy=TRUE, alternative="greater") [2]
}
pop_2_moran <- pop_2_moran[-1,] # Delete empty first line
pop_2_moran_round <- pop_2_moran %>% mutate(across(starts_with("pop"), round, 3))
pop_2_moran_roundCombine Moran tables pop_0, pop_1, pop_2
pop_all_moran <- bind_cols(pop_0_moran_round, pop_1_moran_round, pop_2_moran_round)Analysis of empirical data about dietary habits (section 2.2)
Summary by sociodemographic groups (cf. table 1)
data <- read.csv("data/sociodemo_distribution_healthy.csv", stringsAsFactors = F)
head(data)## Sex Age Edu n_idf N_2002 N_2008 conso_5_2002 conso_5_2008
## 1 1 1 1 493871 52 148 0.01923077 0.02702703
## 2 1 1 2 496964 97 223 0.03092784 0.04932735
## 3 1 1 3 199476 84 137 0.03571429 0.07299270
## 4 1 2 1 1077670 241 355 0.04979253 0.08169014
## 5 1 2 2 662407 122 225 0.07377049 0.11111111
## 6 1 2 3 658465 118 180 0.05932203 0.16666667n_idf = number of individuals of a particular group in
2012 in the Paris Region (Ile de France). N_2002 = number
of individuals of a particular group in the Health and Nutrition
Barometer survey in 2002. N_2008 = number of individuals of
a particular group in the Health and Nutrition Barometer in 2008.
conso_5_2002 = proportion of individuals eating 5+ portions
of fruit and vegetables a day in the Health and Nutrition Barometer in
2002. conso_5_2008 = proportion of individuals eating 5+
portions of fruit and vegetables a day in the Health and Nutrition
Barometer in 2008.
Computation of the total proportion of healthy individuals (i.e. eating 5+ portions of fruit and vegetables a day)
data$conso_5_2002 <- as.numeric(data$conso_5_2002)
data$conso_5_2008 <- as.numeric(data$conso_5_2008)
n_tot_idf <- sum(data$n_idf)
data$share_idf <- data$n_idf / n_tot_idf
data$h_idf_2002 <- round(data$n_idf * data$conso_5_2002,0)
data$h_idf_2008 <- round(data$n_idf * data$conso_5_2008,0)
data$h_baro_2002 <- round(data$N_2002 * data$conso_5_2002,0)
data$h_baro_2008 <- round(data$N_2008 * data$conso_5_2008,0)
total_prop <- function(data, year, sample){
if(sample == "idf"){
tp <-
sum(data[,paste("h", sample, year, sep="_")]) /
sum(data$n_idf)
} else {
tp <-
sum(data[,paste("h", sample, year, sep="_")]) /
sum(data[,paste("N", year, sep="_")])
}
return(tp)
}Results with the Paris region as reference population
total_prop(data=data, year = 2002, sample = "idf")## [1] 0.09568032total_prop(data=data, year = 2008, sample = "idf")## [1] 0.1205141Results with the diet survey as reference population
total_prop(data=data, year = 2002, sample = "baro")## [1] 0.1063931total_prop(data=data, year = 2008, sample = "baro")## [1] 0.1201574Computation of EducIneqIndex
n_tot_baro_2002 <- sum(data$N_2002)
n_tot_baro_2008 <- sum(data$N_2008)
EII <- function(data, year, sample){
if(sample == "idf"){
eii <- (
data[data$Sex == 1 & data$Age == 1 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 1 & data$Age == 1 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 1 & data$Age == 1, "n_idf"] / n_tot_idf)
) + (
data[data$Sex == 1 & data$Age == 2 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 1 & data$Age == 2 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 1 & data$Age == 2, "n_idf"] / n_tot_idf)
) + (
data[data$Sex == 1 & data$Age == 3 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 1 & data$Age == 3 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 1 & data$Age == 3, "n_idf"] / n_tot_idf)
) + (
data[data$Sex == 2 & data$Age == 1 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 2 & data$Age == 1 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 2 & data$Age == 1, "n_idf"] / n_tot_idf)
) + (
data[data$Sex == 2 & data$Age == 2 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 2 & data$Age == 2 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 2 & data$Age == 2, "n_idf"] / n_tot_idf)
) + (
data[data$Sex == 2 & data$Age == 3 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 2 & data$Age == 3 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 2 & data$Age == 3, "n_idf"] / n_tot_idf)
)
} else {
if(year == 2002){ n_tot_baro <- n_tot_baro_2002 }
if(year == 2008){ n_tot_baro <- n_tot_baro_2008 }
eii <- (
data[data$Sex == 1 & data$Age == 1 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 1 & data$Age == 1 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 1 & data$Age == 1, paste0("N_", year)] / n_tot_baro)
) + (
data[data$Sex == 1 & data$Age == 2 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 1 & data$Age == 2 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 1 & data$Age == 2, paste0("N_", year)] / n_tot_baro)
) + (
data[data$Sex == 1 & data$Age == 3 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 1 & data$Age == 3 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 1 & data$Age == 3, paste0("N_", year)] / n_tot_baro)
) + (
data[data$Sex == 2 & data$Age == 1 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 2 & data$Age == 1 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 2 & data$Age == 1, paste0("N_", year)] / n_tot_baro)
) + (
data[data$Sex == 2 & data$Age == 2 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 2 & data$Age == 2 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 2 & data$Age == 2, paste0("N_", year)] / n_tot_baro)
) + (
data[data$Sex == 2 & data$Age == 3 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 2 & data$Age == 3 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 2 & data$Age == 3, paste0("N_", year)] / n_tot_baro)
)
}
return(eii)
}Inequality results with the Paris region as reference population
EII(data=data, year = 2002, sample = "idf")## [1] 1.407512EII(data=data, year = 2008, sample = "idf")## [1] 1.89476Inequality results with the dietary survey as reference population
EII(data=data, year = 2002, sample = "baro")## [1] 1.382566EII(data=data, year = 2008, sample = "baro")## [1] 1.905816Evaluation of 5aday ABM model (section 3.4)
Importing calibration results
pareto <- read.csv2("data/nsga2_clean.csv")
head(pareto)## evolution.generation evolution.evaluated evolution.samples maxProbaToSwitch
## 1 376800 376800 1 0.751
## 2 376800 376800 6 0.750
## 3 376800 376800 2 0.751
## 4 376800 376800 9 0.749
## 5 376800 376800 1 0.749
## 6 376800 376800 1 0.748
## constraintsStrength inertiaCoefficient healthyDietReward
## 1 0.100 1 0
## 2 0.100 1 0
## 3 0.098 1 0
## 4 0.097 1 0
## 5 0.097 1 0
## 6 0.097 1 0
## objective.deltaHealth objective.deltaSocialInequality X.Healthy_vs._2002
## 1 331,736,427 0.138 down
## 2 333,771,427 0.149 down
## 3 334,153,927 0.140 down
## 4 334,658,427 0.128 down
## 5 335,441,427 0.123 down
## 6 335,744,427 0.123 down
## X.Healthy_vs._2008 socialIneq_vs._2002 typo typo_name
## 1 down up down_down_up decrease%Healthy
## 2 down up down_down_up decrease%Healthy
## 3 down up down_down_up decrease%Healthy
## 4 down up down_down_up decrease%Healthy
## 5 down up down_down_up decrease%Healthy
## 6 down up down_down_up decrease%Healthy
## socialInequality
## 1 [1.7569191058549687]
## 2 [1.7525945924779105,1.7540192352458608,1.7427673698247963,1.749312240252607,1.7412857735082057,1.742133022717344]
## 3 [1.7625539632685843,1.745986394640198]
## 4 [1.7662218959212503,1.7696432627743401,1.7665457668386697,1.756434185222282,1.7667110959021646,1.7696451049204427,1.7694793140516023,1.763192441014371,1.7589183006338713]
## 5 [1.7713301959795948]
## 6 [1.7720880068092864]
## deltaSocialInequality
## 1 [0.13782823564503133]
## 2 [0.14215274902208952,0.14072810625413923,0.15197997167520372,0.14543510124739312,0.15346156799179433,0.15261431878265608]
## 3 [0.13219337823141575,0.14876094685980212]
## 4 [0.12852544557874968,0.12510407872565987,0.1282015746613303,0.138313156277718,0.12803624559783544,0.12510223657955732,0.1252680274483977,0.13155490048562912,0.13582904086612868]
## 5 [0.12341714552040517]
## 6 [0.12265933469071366]
## opinionMedian
## 1 [0.49988672137260437]
## 2 [0.4999333918094635,0.4999331533908844,0.49993330240249634,0.49993327260017395,0.4999333322048187,0.49993306398391724]
## 3 [0.4997502863407135,0.49975064396858215]
## 4 [0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5]
## 5 [0.5]
## 6 [0.499999463558197]
## numberOfHealthy
## 1 [834951]
## 2 [834948,832940,835222,834168,833731,834948]
## 3 [827118,826925]
## 4 [818522,817720,816927,815986,817833,817780,815267,817438,818191]
## 5 [817706]
## 6 [811529]
## deltaHealth
## 1 [331736.42711645283]
## 2 [332013.5408605285,334701.42711645283,332841.42711645283,332113.42711645283,335052.42711645283,334725.42711645283]
## 3 [333369.42711645283,334938.42711645283]
## 4 [333959.42711645283,333931.42711645283,335330.42711645283,334215.42711645283,336478.42711645283,336057.42711645283,334658.42711645283,335249.42711645283,334452.42711645283]
## 5 [335441.42711645283]
## 6 [335744.42711645283]Drawing the Pareto front to select a parameter set (figure 4)
pareto$reliability.sample <- ifelse(pareto$evolution.samples>5, 2,1)
pareto$objective.deltaHealth <- as.numeric(gsub(",","",pareto$objective.deltaHealth))
front <- ggplot(data=pareto, aes(x = objective.deltaHealth, y=objective.deltaSocialInequality))
front +
geom_point(aes(group = typo_name, colour = typo_name,
size=reliability.sample)) +
scale_size_continuous(range = c(0.5,1.5)) +
scale_colour_manual(values = alpha(
c("#faa911", #orange
"#07b55e", #green
"#1088b0", #blue
"#aae860" #light green
), .5)) +
geom_vline(xintercept = 290000000) +
geom_hline(yintercept = 0.36) +
theme_bw() 
Distribution of parameter values along the Pareto front (Appendix 6)
pareto$id <- rownames(pareto)
colnames(pareto)## [1] "evolution.generation" "evolution.evaluated"
## [3] "evolution.samples" "maxProbaToSwitch"
## [5] "constraintsStrength" "inertiaCoefficient"
## [7] "healthyDietReward" "objective.deltaHealth"
## [9] "objective.deltaSocialInequality" "X.Healthy_vs._2002"
## [11] "X.Healthy_vs._2008" "socialIneq_vs._2002"
## [13] "typo" "typo_name"
## [15] "socialInequality" "deltaSocialInequality"
## [17] "opinionMedian" "numberOfHealthy"
## [19] "deltaHealth" "reliability.sample"
## [21] "id"pareto_long <- melt(pareto[,c("maxProbaToSwitch","constraintsStrength" ,
"inertiaCoefficient","healthyDietReward",
"objective.deltaSocialInequality", "objective.deltaHealth",
"typo_name", "reliability.sample")],
id.vars=c("objective.deltaSocialInequality", "objective.deltaHealth",
"typo_name", "reliability.sample"))
ggplot(pareto_long, aes(objective.deltaHealth, value)) +
geom_point(aes(group = typo_name, colour = typo_name,
size=reliability.sample)) +
scale_size_continuous(range = c(0.5,1.5)) + scale_colour_manual(values = alpha(
c("#faa911", #orange
"#07b55e", #green
"#1088b0", #blue
"#aae860" #light green
), .5)) +
theme_bw() + theme(legend.position="none") +
facet_wrap(~variable, scales = "free_y", ncol = 1)
ggplot(pareto_long, aes(value, objective.deltaSocialInequality)) +
geom_point(aes(group = typo_name, colour = typo_name,
size=reliability.sample)) +
scale_size_continuous(range = c(0.5,1.5)) + scale_colour_manual(values = alpha(
c("#faa911", #orange
"#07b55e", #green
"#1088b0", #blue
"#aae860" #light green
), .5)) +
theme_bw() + theme(legend.position="none")+
facet_wrap(~variable, scales = "free_x", nrow = 1)
Results (section 4)
Spatial distribution of agents with healthy behaviour for the 5 scenarios - in the night cells (section 4.3)
s0r <- st_read( "data/results_HigherProp_Scenario1_RandomPop_NoMove_42.gpkg",
layer="Initial_0", quiet = TRUE)
s0o <- st_read( "data/results_HigherProp_Scenario5_ObservedPop_ObservedMove_42.gpkg",
layer="Initial_0", quiet = TRUE)
s1A <- st_read( "data/results_HigherProp_Scenario1_RandomPop_NoMove_42.gpkg",
layer="Final", quiet = TRUE)
s1B <- st_read( "data/results_HigherProp_Scenario2_RandomPop_RandomMove_42.gpkg",
layer="Final", quiet = TRUE)
s2A <- st_read( "data/results_HigherProp_Scenario3_ObservedPop_NoMove_42.gpkg",
layer="Final", quiet = TRUE)
s2B <- st_read( "data/results_HigherProp_Scenario4_ObservedPop_RandomMove_42.gpkg",
layer="Final", quiet = TRUE)
s2C <- st_read( "data/results_HigherProp_Scenario5_ObservedPop_ObservedMove_42.gpkg",
layer="Final", quiet = TRUE)
bks <- c(0, 0.1, 0.125, 0.15, 0.175, Inf)Maps (figure 6)
Initial Random
tm_shape(s0r) + tm_fill(col="X.propHealthy", style="fixed",
breaks = bks, palette="Greens") +
tm_legend(outside=TRUE)
Initial Observed
tm_shape(s0o) + tm_fill(col="X.propHealthy", style="fixed",
breaks = bks, palette="Greens") +
tm_legend(outside=TRUE)
Scenario 1A
tm_shape(s1A) + tm_fill(col="X.propHealthy", style="fixed",
breaks = bks, palette="Greens") +
tm_legend(outside=TRUE)
Scenario 1B
tm_shape(s1B) + tm_fill(col="X.propHealthy", style="fixed",
breaks = bks, palette="Greens") +
tm_legend(outside=TRUE)
Scenario 2A
tm_shape(s2A) + tm_fill(col="X.propHealthy", style="fixed",
breaks = bks, palette="Greens") +
tm_legend(outside=TRUE)
Scenario 2B
tm_shape(s2B) + tm_fill(col="X.propHealthy", style="fixed",
breaks = bks, palette="Greens") +
tm_legend(outside=TRUE)
Scenario 2C
tm_shape(s2C) + tm_fill(col="X.propHealthy", style="fixed",
breaks = bks, palette="Greens") +
tm_legend(outside=TRUE)
Duncan Dissimilarity Index (figure 6)
duncan <- data.frame(matrix(ncol = 2, nrow=7))
colnames(duncan) <- c("type", "Duncan" )
dfs0r<- s0r %>% st_drop_geometry()
duncan[1,1] <- "s0r_night"
duncan[1,2]<- ISDuncan(dfs0r[,3:4]) [1]
dfs0o<- s0o %>% st_drop_geometry()
duncan[2,1] <- "s0o_night"
duncan[2,2]<- ISDuncan(dfs0o[,3:4]) [1]
dfs1A<- s1A %>% st_drop_geometry()
duncan[3,1] <- "s1A_final_night"
duncan[3,2]<- ISDuncan(dfs1A[,3:4]) [1]
dfs1B<- s1B %>% st_drop_geometry()
duncan[4,1] <- "s1B_final_night"
duncan[4,2]<- ISDuncan(dfs1B[,3:4]) [1]
dfs2A<- s2A %>% st_drop_geometry()
duncan[5,1] <- "s2A_final_night"
duncan[5,2]<- ISDuncan(dfs2A[,3:4]) [1]
dfs2B<- s2B %>% st_drop_geometry()
duncan[6,1] <- "s2B_final_night"
duncan[6,2]<- ISDuncan(dfs2B[,3:4]) [1]
dfs2C<- s2C %>% st_drop_geometry()
duncan[7,1] <- "s2C_final_night"
duncan[7,2]<- ISDuncan(dfs2C[,3:4]) [1]
duncan_round <- duncan %>% mutate(across(starts_with("Duncan"), round, 3))
duncan_roundMoran’s index of spatial autocorrelation (figure 6)
moran <- data.frame(matrix(ncol = 3, nrow=7))
colnames(moran) <- c("type", "IMoran" , "pvalueMoran")
nbs0r <- poly2nb(s0r, queen=TRUE)
lws0r <- nb2listw(nbs0r, style="W", zero.policy=TRUE)
moran[1,1] <- "s0r_night"
moran[1,2] <- moran.test(s0r$X.propHealthy, lws0r, zero.policy=TRUE, alternative="greater") [3]
moran[1,3] <- moran.test(s0r$X.propHealthy, lws0r, zero.policy=TRUE, alternative="greater") [2]
nbs0o <- poly2nb(s0o, queen=TRUE)
lws0o <- nb2listw(nbs0o, style="W", zero.policy=TRUE)
moran[4,1] <- "s0o_night"
moran[4,2] <- moran.test(s0o$X.propHealthy, lws0o, zero.policy=TRUE, alternative="greater") [3]
moran[4,3] <- moran.test(s0o$X.propHealthy, lws0o, zero.policy=TRUE, alternative="greater") [2]
nbs1A <- poly2nb(s1A, queen=TRUE)
lws1A <- nb2listw(nbs1A, style="W", zero.policy=TRUE)
moran[2,1] <- "s1A_final_night"
moran[2,2] <- moran.test(s1A$X.propHealthy, lws1A, zero.policy=TRUE, alternative="greater") [3]
moran[2,3] <- moran.test(s1A$X.propHealthy, lws1A, zero.policy=TRUE, alternative="greater") [2]
nbs1B <- poly2nb(s1B, queen=TRUE)
lws1B <- nb2listw(nbs1B, style="W", zero.policy=TRUE)
moran[3,1] <- "s1B_final_night"
moran[3,2] <- moran.test(s1B$X.propHealthy, lws1B, zero.policy=TRUE, alternative="greater") [3]
moran[3,3] <- moran.test(s1B$X.propHealthy, lws1B, zero.policy=TRUE, alternative="greater") [2]
nbs2A <- poly2nb(s2A, queen=TRUE)
lws2A <- nb2listw(nbs2A, style="W", zero.policy=TRUE)
moran[5,1] <- "s2A_final_night"
moran[5,2] <- moran.test(s2A$X.propHealthy, lws2A, zero.policy=TRUE, alternative="greater") [3]
moran[5,3] <- moran.test(s2A$X.propHealthy, lws2A, zero.policy=TRUE, alternative="greater") [2]
nbs2B <- poly2nb(s2B, queen=TRUE)
lws2B <- nb2listw(nbs2B, style="W", zero.policy=TRUE)
moran[6,1] <- "s2B_final_night"
moran[6,2] <- moran.test(s2B$X.propHealthy, lws2B, zero.policy=TRUE, alternative="greater") [3]
moran[6,3] <- moran.test(s2B$X.propHealthy, lws2B, zero.policy=TRUE, alternative="greater") [2]
nbs2C <- poly2nb(s2C, queen=TRUE)
lws2C <- nb2listw(nbs2C, style="W", zero.policy=TRUE)
moran[7,1] <- "s2C_final_night"
moran[7,2] <- moran.test(s2C$X.propHealthy, lws2C, zero.policy=TRUE, alternative="greater") [3]
moran[7,3] <- moran.test(s2C$X.propHealthy, lws2C, zero.policy=TRUE, alternative="greater") [2]
moran_round <- moran %>% mutate(across(2:3, round, 3))
moran_round ## type IMoran pvalueMoran
## 1 s0r_night 0.010 0.061
## 2 s1A_final_night 0.010 0.057
## 3 s1B_final_night 0.016 0.008
## 4 s0o_night 0.006 0.175
## 5 s2A_final_night 0.033 0.000
## 6 s2B_final_night 0.024 0.000
## 7 s2C_final_night 0.009 0.091Alternative measure of social inequality in health behaviour (Appendix 4)
Computation of rank-ordered inequality measure
data_edu <- aggregate(data[,c("n_idf", "N_2002", "N_2008",
"h_idf_2002", "h_idf_2008",
"h_baro_2002", "h_baro_2008")],
by = list(data$Edu), FUN = sum)
colnames(data_edu)[1] <- "edu"
data_edu## edu n_idf N_2002 N_2008 h_idf_2002 h_idf_2008 h_baro_2002 h_baro_2008
## 1 1 3800999 1030 1372 336079 359509 122 155
## 2 2 2946028 595 1135 279323 361858 57 132
## 3 3 2021871 471 797 223609 335409 44 110guido <- function(data_edu, year, sample){
if(sample == "idf"){
n <- sum(data_edu[,"n_idf"])
n1 <- data_edu[data_edu$edu == 1,"n_idf"]
n2 <- data_edu[data_edu$edu == 2,"n_idf"]
n3 <- data_edu[data_edu$edu == 3,"n_idf"]
} else {
n <- sum(data_edu[,paste("N", year, sep="_")])
n1 <- data_edu[data_edu$edu == 1,paste("N", year, sep="_")]
n2 <- data_edu[data_edu$edu == 2,paste("N", year, sep="_")]
n3 <- data_edu[data_edu$edu == 3,paste("N", year, sep="_")]
}
h1 <- data_edu[data_edu$edu == 1,paste("h", sample, year, sep="_")]
h2 <- data_edu[data_edu$edu == 2,paste("h", sample, year, sep="_")]
h3 <- data_edu[data_edu$edu == 3,paste("h", sample, year, sep="_")]
t1 <- (1 - n2 - n3) / 2
t2 <- (n1 - n3 + 1) / 2
t3 <- (n1 + n2 + 1) / 2
guido <- ( 8 / n^2 ) * (h1 * t1 + h2 * t2 + h3 * t3)
return(guido)
}
guido(data_edu = data_edu, year = 2002, sample = "idf")## [1] 0.01748087guido(data_edu = data_edu, year = 2008, sample = "idf")## [1] 0.05830406guido(data_edu = data_edu, year = 2002, sample = "baro")## [1] -0.02409715guido(data_edu = data_edu, year = 2008, sample = "baro")## [1] 0.01927629Density distribution of rank-order inequality index per scenario
mean_1A <- mean(sim[sim$numOfScenario == 1,"erreygersE"])
mean_1B <- mean(sim[sim$numOfScenario == 2,"erreygersE"])
mean_2A <- mean(sim[sim$numOfScenario == 3,"erreygersE"])
mean_2B <- mean(sim[sim$numOfScenario == 4,"erreygersE"])
mean_2C <- mean(sim[sim$numOfScenario == 5,"erreygersE"])
sim %>%
ggplot( aes(x=erreygersE, group=Scenario, fill=Scenario)) +
geom_vline(xintercept=0.0175, size=1, color="grey30", linetype = "dashed") +
geom_density(alpha=0.6) +
scale_fill_manual(values=c("firebrick2", "mediumorchid4", "dodgerblue2", "darkorange2", "springgreen3"))+
theme_bw() +
geom_vline(xintercept=mean_1A, size=0.4, color="firebrick2") +
geom_vline(xintercept=mean_1B, size=0.4, color="mediumorchid4") +
geom_vline(xintercept=mean_2A, size=0.4, color="dodgerblue2") +
geom_vline(xintercept=mean_2B, size=0.4, color="darkorange2") +
geom_vline(xintercept=mean_2C, size=0.4, color="springgreen3") +
annotate("text", x=mean_1A - 0.004, y= -10, size=2.5, label=round(mean_1A,3), color="firebrick2") +
annotate("text", x=mean_1B - 0.004 , y= -10, size=2.5, label=round(mean_1B,3), color="mediumorchid4") +
annotate("text", x=mean_2A - 0.004 , y= -10, size=2.5, label=round(mean_2A,3), color="dodgerblue2") +
annotate("text", x=mean_2B + 0.004 , y= -10, size=2.5, label=round(mean_2B,3), color="darkorange2") +
annotate("text", x=mean_2C + 0.004, y= -10, size=2.5, label=round(mean_2C,3), color="springgreen3") +
annotate("text", x=0.0175 + 0.0051, y= 115, size=3, label="Value at initialisation", color="grey30") +
xlab("Rank-Order Inequality Index") +
ylab("Number of simulations") +
guides(fill=guide_legend(title="Type of scenario")) #+ 
Alternative measures of spatial distribution of agents with healthy behaviour - in the day and evening cells (Appendix 8)
In day cells
s1A_day <- st_read( "data/results_HigherProp_Scenario1_RandomPop_NoMove_42.gpkg", quiet = TRUE,
layer="Final_1")
s1B_day <- st_read( "data/results_HigherProp_Scenario2_RandomPop_RandomMove_42.gpkg", quiet = TRUE,
layer="Final_1")
s2A_day <- st_read( "data/results_HigherProp_Scenario3_ObservedPop_NoMove_42.gpkg", quiet = TRUE,
layer="Final_1")
s2B_day <- st_read( "data/results_HigherProp_Scenario4_ObservedPop_RandomMove_42.gpkg", quiet = TRUE,
layer="Final_1")
s2C_day <- st_read( "data/results_HigherProp_Scenario5_ObservedPop_ObservedMove_42.gpkg", quiet = TRUE,
layer="Final_1")Duncan
duncan_day <- data.frame(matrix(ncol = 2, nrow=5))
colnames(duncan_day) <- c("type", "Duncan_dayCell" )
dfs1A_day<- s1A_day %>% st_drop_geometry()
duncan_day[1,1] <- "s1A_final_day"
duncan_day[1,2]<- ISDuncan(dfs1A_day[,3:4]) [1]
dfs1B_day<- s1B_day %>% st_drop_geometry()
duncan_day[2,1] <- "s1B_final_day"
duncan_day[2,2]<- ISDuncan(dfs1B_day[,3:4]) [1]
dfs2A_day<- s2A_day %>% st_drop_geometry()
duncan_day[3,1] <- "s2A_final_day"
duncan_day[3,2]<- ISDuncan(dfs2A_day[,3:4]) [1]
dfs2B_day<- s2B_day %>% st_drop_geometry()
duncan_day[4,1] <- "s2B_final_day"
duncan_day[4,2]<- ISDuncan(dfs2B_day[,3:4]) [1]
dfs2C_day<- s2C_day %>% st_drop_geometry()
duncan_day[5,1] <- "s2C_final_day"
duncan_day[5,2]<- ISDuncan(dfs2C_day[,3:4]) [1]
duncan_day <- duncan_day %>% mutate(across(starts_with("Duncan"), round, 3))Moran
moran_day <- data.frame(matrix(ncol = 3, nrow=5))
colnames(moran_day) <- c("type", "IMoran_dayCell" , "pvalueMoran_dayCell")
nbs1A_day <- poly2nb(s1A_day, queen=TRUE)
lws1A_day <- nb2listw(nbs1A_day, style="W", zero.policy=TRUE)
moran_day[1,1] <- "s1A_final_day"
moran_day[1,2] <- moran.test(s1A_day$X.propHealthy, lws1A_day, zero.policy=TRUE, alternative="greater") [3]
moran_day[1,3] <- moran.test(s1A_day$X.propHealthy, lws1A_day, zero.policy=TRUE, alternative="greater") [2]
nbs1B_day <- poly2nb(s1B_day, queen=TRUE)
lws1B_day <- nb2listw(nbs1B_day, style="W", zero.policy=TRUE)
moran_day[2,1] <- "s1B_final_day"
moran_day[2,2] <- moran.test(s1B_day$X.propHealthy, lws1B_day, zero.policy=TRUE, alternative="greater") [3]
moran_day[2,3] <- moran.test(s1B_day$X.propHealthy, lws1B_day, zero.policy=TRUE, alternative="greater") [2]
nbs2A_day <- poly2nb(s2A_day, queen=TRUE)
lws2A_day <- nb2listw(nbs2A_day, style="W", zero.policy=TRUE)
moran_day[3,1] <- "s2A_final_day"
moran_day[3,2] <- moran.test(s2A_day$X.propHealthy, lws2A_day, zero.policy=TRUE, alternative="greater") [3]
moran_day[3,3] <- moran.test(s2A_day$X.propHealthy, lws2A_day, zero.policy=TRUE, alternative="greater") [2]
nbs2B_day <- poly2nb(s2B_day, queen=TRUE)
lws2B_day <- nb2listw(nbs2B_day, style="W", zero.policy=TRUE)
moran_day[4,1] <- "s2B_final_day"
moran_day[4,2] <- moran.test(s2B_day$X.propHealthy, lws2B_day, zero.policy=TRUE, alternative="greater") [3]
moran_day[4,3] <- moran.test(s2B_day$X.propHealthy, lws2B_day, zero.policy=TRUE, alternative="greater") [2]
nbs2C_day <- poly2nb(s2C_day, queen=TRUE)
lws2C_day <- nb2listw(nbs2C_day, style="W", zero.policy=TRUE)
moran_day[5,1] <- "s2C_final_day"
moran_day[5,2] <- moran.test(s2C_day$X.propHealthy, lws2C_day, zero.policy=TRUE, alternative="greater") [3]
moran_day[5,3] <- moran.test(s2C_day$X.propHealthy, lws2C_day, zero.policy=TRUE, alternative="greater") [2]
moran_round_day <- moran_day %>% mutate(across(2:3, round, 3))Combine Duncan and Moran files
Duncan_moran_dayCell <- left_join(duncan_day, moran_round_day, by = "type")
Duncan_moran_dayCell In evening cells
s1A_evening <- st_read( "data/results_HigherProp_Scenario1_RandomPop_NoMove_42.gpkg", quiet = TRUE,
layer="Final_2")
s1B_evening <- st_read( "data/results_HigherProp_Scenario2_RandomPop_RandomMove_42.gpkg", quiet = TRUE,
layer="Final_2")
s2A_evening <- st_read( "data/results_HigherProp_Scenario3_ObservedPop_NoMove_42.gpkg", quiet = TRUE,
layer="Final_2")
s2B_evening <- st_read( "data/results_HigherProp_Scenario4_ObservedPop_RandomMove_42.gpkg", quiet = TRUE,
layer="Final_2")
s2C_evening <- st_read( "data/results_HigherProp_Scenario5_ObservedPop_ObservedMove_42.gpkg", quiet = TRUE,
layer="Final_2")Duncan
duncan_evening <- data.frame(matrix(ncol = 2, nrow=5))
colnames(duncan_evening) <- c("type", "Duncan_eveningCell" )
dfs1A_evening<- s1A_evening %>% st_drop_geometry()
duncan_evening[1,1] <- "s1A_final_evening"
duncan_evening[1,2]<- ISDuncan(dfs1A_evening[,3:4]) [1]
dfs1B_evening<- s1B_evening %>% st_drop_geometry()
duncan_evening[2,1] <- "s1B_final_evening"
duncan_evening[2,2]<- ISDuncan(dfs1B_evening[,3:4]) [1]
dfs2A_evening<- s2A_evening %>% st_drop_geometry()
duncan_evening[3,1] <- "s2A_final_evening"
duncan_evening[3,2]<- ISDuncan(dfs2A_evening[,3:4]) [1]
dfs2B_evening<- s2B_evening %>% st_drop_geometry()
duncan_evening[4,1] <- "s2B_final_evening"
duncan_evening[4,2]<- ISDuncan(dfs2B_evening[,3:4]) [1]
dfs2C_evening<- s2C_evening %>% st_drop_geometry()
duncan_evening[5,1] <- "s2C_final_evening"
duncan_evening[5,2]<- ISDuncan(dfs2C_evening[,3:4]) [1]
duncan_evening <- duncan_evening %>% mutate(across(starts_with("Duncan"), round, 3))Moran
moran_evening <- data.frame(matrix(ncol = 3, nrow=5))
colnames(moran_evening) <- c("type", "IMoran_eveningCell" , "pvalueMoran_eveningCell")
nbs1A_evening <- poly2nb(s1A_evening, queen=TRUE)
lws1A_evening <- nb2listw(nbs1A_evening, style="W", zero.policy=TRUE)
moran_evening[1,1] <- "s1A_final_evening"
moran_evening[1,2] <- moran.test(s1A_evening$X.propHealthy, lws1A_evening, zero.policy=TRUE, alternative="greater") [3]
moran_evening[1,3] <- moran.test(s1A_evening$X.propHealthy, lws1A_evening, zero.policy=TRUE, alternative="greater") [2]
nbs1B_evening <- poly2nb(s1B_evening, queen=TRUE)
lws1B_evening <- nb2listw(nbs1B_evening, style="W", zero.policy=TRUE)
moran_evening[2,1] <- "s1B_final_evening"
moran_evening[2,2] <- moran.test(s1B_evening$X.propHealthy, lws1B_evening, zero.policy=TRUE, alternative="greater") [3]
moran_evening[2,3] <- moran.test(s1B_evening$X.propHealthy, lws1B_evening, zero.policy=TRUE, alternative="greater") [2]
nbs2A_evening <- poly2nb(s2A_evening, queen=TRUE)
lws2A_evening <- nb2listw(nbs2A_evening, style="W", zero.policy=TRUE)
moran_evening[3,1] <- "s2A_final_evening"
moran_evening[3,2] <- moran.test(s2A_evening$X.propHealthy, lws2A_evening, zero.policy=TRUE, alternative="greater") [3]
moran_evening[3,3] <- moran.test(s2A_evening$X.propHealthy, lws2A_evening, zero.policy=TRUE, alternative="greater") [2]
nbs2B_evening <- poly2nb(s2B_evening, queen=TRUE)
lws2B_evening <- nb2listw(nbs2B_evening, style="W", zero.policy=TRUE)
moran_evening[4,1] <- "s2B_final_evening"
moran_evening[4,2] <- moran.test(s2B_evening$X.propHealthy, lws2B_evening, zero.policy=TRUE, alternative="greater") [3]
moran_evening[4,3] <- moran.test(s2B_evening$X.propHealthy, lws2B_evening, zero.policy=TRUE, alternative="greater") [2]
nbs2C_evening <- poly2nb(s2C_evening, queen=TRUE)
lws2C_evening <- nb2listw(nbs2C_evening, style="W", zero.policy=TRUE)
moran_evening[5,1] <- "s2C_final_evening"
moran_evening[5,2] <- moran.test(s2C_evening$X.propHealthy, lws2C_evening, zero.policy=TRUE, alternative="greater") [3]
moran_evening[5,3] <- moran.test(s2C_evening$X.propHealthy, lws2C_evening, zero.policy=TRUE, alternative="greater") [2]
moran_round_evening <- moran_evening %>% mutate(across(2:3, round, 3))Combine Duncan and Moran files
Duncan_moran_eveningCell <- left_join(duncan_evening, moran_round_evening, by = "type")
Duncan_moran_eveningCellEvolution over the six series of three time slices in simulated values - 5aDay model (Appendix 9)
library(cowplot)
library(tidyverse)
library(here)phase <- read.csv("/Users/ccottineau/GitHub/Paper/Notebook_R_fig_tables/data/phase_48.csv")
long_table <- data.frame()
i <- 0
j <- 1
time <- data.frame(
day = rep("day",18),
nd = sort(rep(0:5,3)),
step = rep("step",18),
ns = rep(0:2,6)) |>
mutate(time = paste0(day, nd, "_", step, ns)) |>
select(time)
for(i in 1:nrow(phase)){
nsenario <- phase[i, "numOfScenario"]
EII <- t(
str_split(
str_replace(
str_replace( phase[i, "socialInequality"] , c("\\["), ""),
c("\\["), ""
),
",", simplify = T
)
)[1:18,]
healthy <- t(
str_split(
str_replace(
str_replace( phase[i, "numberOfHealthy"] , c("\\["), ""),
c("\\["), ""
),
",", simplify = T
)
)[1:18,]
long_table[j:(j+17),"time"] <- time
long_table[j:(j+17),"scenario"] <- rep(nsenario, 18)
long_table[j:(j+17),"socialInequality"] <- EII
long_table[j:(j+17),"numberOfHealthy"] <- healthy
j <- j+18
}
plot_table <- long_table |>
mutate(prop_healthy = as.numeric(numberOfHealthy) / 8768898,
EII = as.numeric(socialInequality),
scenario = as.factor(ifelse(scenario == 2, "1B",
ifelse(scenario == 5, "2C", NA))
)
)
ggplot(plot_table,
aes(x = time, y = EII, group = scenario,
colour = scenario)) +
geom_point()
specialRead <- function(x) {
read_csv(file=x, col_names = TRUE, col_types= "nnnnnnnnnn")
}
concatFiles <- function (base, folder){
fullPath = paste(base,folder,sep="")
list.files(path = fullPath, full.names = TRUE, pattern = "(\\d_\\d_\\d)") %>% lapply(specialRead) %>% bind_rows
}
recode <- function(x) {
step_to_day <- c (
'00' = 'day 0 step 1',
'01' = 'day 0 step 2',
'02' = 'day 0 step 3',
'10' = 'day 1 step 1',
'11' = 'day 1 step 2',
'12' = 'day 1 step 3',
'20' = 'day 2 step 1',
'21' = 'day 2 step 2',
'22' = 'day 2 step 3',
'30' = 'day 3 step 1',
'31' = 'day 3 step 2',
'32' = 'day 3 step 3',
'40' = 'day 4 step 1',
'41' = 'day 4 step 2',
'42' = 'day 4 step 3',
'50' = 'day 5 step 1',
'51' = 'day 5 step 2',
'52' = 'day 5 step 3'
)
return(step_to_day[x])
}Simulated proportions of healthy agents, of social inequality (EII) and of rank-ordered inequalities (Et)
hp_sc5_health <- read_csv("data/HigherProp/health.csv")
hp_sc5_health$propHealthy <- hp_sc5_health$healthy / hp_sc5_health$effective
hpLong_sc5_health <- hp_sc5_health %>%
gather(key=facteur, value=valeur, propHealthy,avgOpinion,socialInequality,e ) %>%
mutate(ds= recode(paste0(day, slice)))
hpLong_sc5_health <- hpLong_sc5_health %>% mutate (type="hp")
both_sc5_health <- hpLong_sc5_health %>% mutate(type = as.factor(type)) p_hp_propHealthy <- both_sc5_health %>% filter(facteur=="propHealthy") %>% ggplot( aes(x=ds, y=valeur,colour=type) )+
theme_bw()+
theme(axis.text.x = element_text(angle = 45, vjust = 0.9, hjust=1.0)) +
theme(axis.text.x.bottom = element_text(margin=margin(t=0, r=0,b=10,l=0)) ) +
geom_point() +
theme(legend.position = "none") +
labs(y= "Value", x = "prop. of healthy")
p_hp_ineq <- both_sc5_health %>% filter(facteur=="socialInequality") %>% ggplot( aes(x=ds, y=valeur,colour=type)) +
theme_bw()+
theme(axis.text.x = element_text(angle = 45, vjust = 0.9, hjust=1.0)) +
theme(axis.text.x.bottom = element_text(margin=margin(t=0, r=0,b=10,l=0)) ) +
geom_point() +
theme(legend.position = "none") +
labs(y= "Value", x = "EducIneqIndex (EII)")
p_hp_e <- both_sc5_health %>% filter(facteur=="e") %>% ggplot( aes(x=ds, y=valeur,colour=type)) +
theme_bw()+
theme(axis.text.x = element_text(angle = 45, vjust = 0.9, hjust=1.0)) +
theme(axis.text.x.bottom = element_text(margin=margin(t=0, r=0,b=10,l=0)) ) +
geom_point() +
theme(legend.position = "none") +
labs(y= "Value", x = "Rank-ordered inequalities (Et)")
p_hp_all <- plot_grid(p_hp_propHealthy, p_hp_ineq, p_hp_e, labels = "AUTO",nrow=1, align="h")
p_hp_all
##ggsave("out/p_hp_all.pdf",p_hp_all, width=20)Simulated proportion of healthy agents and simulated opinion across the 18 sociodemographic categories
sortie_sc5_higherprop <- concatFiles("data", "/HigherProp")
sortie_sc5_higherprop$propHealthy <- sortie_sc5_higherprop$healthy / sortie_sc5_higherprop$effective
sortieLong_sc5_hp <- gather(sortie_sc5_higherprop, facteur, valeur, c(propHealthy, avgOpinion) )
sortieLong_sc5_hp <- mutate(sortieLong_sc5_hp, "ds" = gsub(" ", "", paste(paste("[",paste(sortieLong_sc5_hp$day, sortieLong_sc5_hp$slice)),"]")))
sortieLong_sc5_hp <- mutate(sortieLong_sc5_hp, ds= recode(paste0(day, slice)))
age.l <- c("15-29 yrs", "30-59 yrs", "60-75 yrs")
names(age.l) <- c("1", "2", "3")
sex.l <- c("Male", "Female")
names(sex.l) <- c("1", "2")
p1 <- ggplot(data=sortieLong_sc5_hp %>% filter(facteur=="propHealthy"), aes(x=ds, y=valeur, group=educ, fill=educ, colour=as.character(educ)))
p1 <- p1 + scale_color_discrete(
name = "Education level",
labels = c("low", "middle", "high")
)
p1 <- p1 + theme_bw() + theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
p1 <- p1 + labs(title="Proportion of healthy by gender, age and education level", y = "Value", x = "Time Slice")
p1 <- p1 + geom_line(size=0.5) + facet_grid( sex ~ age, labeller = labeller(sex = sex.l, age = age.l), scales = "fixed") #show.legend = FALSE
p1 <- p1 + coord_cartesian( ylim = c(0, 0.3))
p2 <- ggplot(data=sortieLong_sc5_hp %>% filter(facteur=="avgOpinion"), aes(x=ds, y=valeur, group=educ, fill=educ, colour=as.character(educ)))
p2 <- p2 + scale_color_discrete(
name = "Education level",
labels = c("low", "middle", "high")
)
p2 <- p2 + theme_bw() + theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
p2 <- p2 + labs(title="Average Opinion by gender, age and education level", color = "Education level", y = "Value", x = "Time Slice")
p2 <- p2 + geom_line(size=0.5) + facet_grid(sex ~ age, labeller = labeller(sex = sex.l, age = age.l), scales = "fixed") #show.legend = FALSE
p2 <- p2 + coord_cartesian( ylim = c(0.3, 0.7))
propHealtyAndOpinion_higherprop <- plot_grid(
p1, p2,
labels = "AUTO", ncol = 1)
print(propHealtyAndOpinion_higherprop)
##ggsave("out/propHealtyAndOpinion_higherprop.pdf")
Social inequalities in healthy behaviour for the 5 scenarios (section 4.1)
Importing simulation results
Distribution of simulated values of EducIneqIndex (EII) for the five space-time scenarios (10 000 replications per scenario) (figure 5)
Overall proportions of agents with healthy behaviour for the 5 scenarios (section 4.2)
!! TO DO JULIE (Table 5)